Overview

IAT trial type D scores are calculated from an average of only 18 pairs of reaction times. This would be deemed as far too low anywhere else in the reaction time literature. The implications of this can be seen in how poorly estimated any one IAT D score is. We can observe this by bootstrapping reaction times for each participant and trial type.

# dependencies
library(tidyverse)
library(knitr)
library(kableExtra)
library(rsample)
library(broom)
library(purrr)
library(furrr)

# function to round all numeric vars in a data frame
round_df <- function(df, n_digits = 3) {
  df %>% mutate_if(is.numeric, round, digits = n_digits)
}

# run furrr:::future_map in parallel
future::plan(multiprocess)
options(future.globals.maxSize = 3000 * 1024^2)

# options
options(knitr.table.format = "html") # necessary configuration of tables

# disable scientific notation
options(scipen = 999) 

# get data 
data_iat_for_scoring_subset <- read_rds("../data/data_iat_for_scoring_subset.rds")

Descriptives

# data_descriptives <- data_for_analysis %>%
#   distinct(session_id, .keep_all = TRUE) 
# 
# data_descriptives %>%
#   count(domain) %>% 
#   kable() %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
# 
# data_descriptives %>%
#   count(domain) %>% 
#   summarize(total_n = sum(n),
#             min_n_per_domain = min(n),
#             max_n_per_domain = max(n),
#             mean_n_per_domain = round(mean(n, na.rm = TRUE), 2),
#             sd_n_per_domain = round(sd(n, na.rm = TRUE), 2)) %>% 
#   kable() %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
# 
# data_descriptives %>%
#   summarize(min_age = round(min(age, na.rm = TRUE), 2),
#             max_age = round(max(age, na.rm = TRUE), 2),
#             mean_age = round(mean(age, na.rm = TRUE), 2),
#             sd_age = round(sd(age, na.rm = TRUE), 2)) %>% 
#   kable() %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
# 
# data_descriptives %>%
#   count(gender) %>% 
#   kable() %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Bootstrap 95% CIs on D scores

Bootstrapping CIs

circa 5 hours runtime

# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_estimates_with_CIs.rds")) {
  
  data_estimates_with_CIs <- read_rds("models/data_estimates_with_CIs.rds")
  
} else {
  
  # n boots for all metrics
  nboots <- 2000 
  
  # trim RTs>10000 ms
  data_trimmed <- data_iat_for_scoring_subset %>%
    filter(trial_latency <= 10000) 
  
  # create D scores
  data_D_scores <- data_trimmed %>%
    group_by(session_id) %>%
    summarise(Da = (mean(trial_latency[block_number == 5]) - mean(trial_latency[block_number == 2])) /
                sd(trial_latency[block_number %in% c(2, 5)]),
              Db = (mean(trial_latency[block_number == 6]) - mean(trial_latency[block_number == 3])) /
                sd(trial_latency[block_number %in% c(3, 6)])) %>% 
    mutate(D = (Da + Db)/2) %>%
    ungroup() %>%
    select(-Da, -Db) %>%
    round_df(3)
  
  # create cohens d scores
  data_d_scores <- data_trimmed %>%
    group_by(session_id) %>%
    summarise(da = (mean(trial_latency[block_number == 5]) - mean(trial_latency[block_number == 2])) /
                sqrt((sd(trial_latency[block_number == 5])^2 + sd(trial_latency[block_number == 2])^2)/2),
              db = (mean(trial_latency[block_number == 6]) - mean(trial_latency[block_number == 3])) /
                sqrt((sd(trial_latency[block_number == 6])^2 + sd(trial_latency[block_number == 3])^2)/2)) %>%
    mutate(d = (da + db)/2) %>%
    ungroup() %>%
    select(-da, -db) %>%
    round_df(3)

  # function to apply to each resample
  calc_D <- function(split) {
    analysis(split) %>%
      group_by(session_id) %>%
      summarise(Da = (mean(trial_latency[block_number == 5]) - mean(trial_latency[block_number == 2])) /
                  sd(trial_latency[block_number %in% c(2, 5)]),
                Db = (mean(trial_latency[block_number == 6]) - mean(trial_latency[block_number == 3])) /
                  sd(trial_latency[block_number %in% c(3, 6)])) %>% 
      mutate(D = (Da + Db)/2) %>%
      ungroup() %>%
      select(-Da, -Db) %>%
      round_df(3)
  }
  
  # function to apply to each resample
  calc_cohens_d <- function(split) {
    analysis(split) %>%
      group_by(session_id) %>%
      summarise(da = (mean(trial_latency[block_number == 5]) - mean(trial_latency[block_number == 2])) /
                  sqrt((sd(trial_latency[block_number == 5])^2 + sd(trial_latency[block_number == 2])^2)/2),
                db = (mean(trial_latency[block_number == 6]) - mean(trial_latency[block_number == 3])) /
                  sqrt((sd(trial_latency[block_number == 6])^2 + sd(trial_latency[block_number == 3])^2)/2)) %>%
      mutate(d = (da + db)/2) %>%
      ungroup() %>%
      select(-da, -db) %>%
      round_df(3)
  }
  
  # start timer
  start <- Sys.time()

  domains <- data_trimmed %>%
    distinct(domain) %>%
    pull(domain)
  
  # apply to each bootstrap
  data_D_bootstrapped_CIs <- 
    dplyr::bind_rows(
      lapply(seq_along(domains), function(i) {
        data_trimmed %>%
          filter(domain == domains[i]) %>%
          group_by(session_id) %>%
          bootstraps(times = nboots) %>%
          mutate(D_metrics = furrr::future_map(splits, calc_D)) %>%
          select(-splits) %>%
          unnest(D_metrics) %>%
          group_by(session_id) %>%
          dplyr::summarize(D_ci_lower = quantile(D, 0.025, na.rm = TRUE),
                           D_ci_upper = quantile(D, 0.975, na.rm = TRUE)) 
      })
    ) %>%
    mutate(D_sig = ifelse((D_ci_lower < 0 & D_ci_upper < 0) | (D_ci_lower > 0 & D_ci_upper > 0), TRUE, FALSE),
           D_ci_width = D_ci_upper - D_ci_lower) %>%
    round_df(3)
  
  data_d_bootstrapped_CIs <- 
    dplyr::bind_rows(
      lapply(seq_along(domains), function(i) {
        data_trimmed %>%
          filter(domain == domains[i]) %>%
          group_by(session_id) %>%
          bootstraps(times = nboots) %>%
          mutate(d_metrics = furrr::future_map(splits, calc_cohens_d)) %>%
          select(-splits) %>%
          unnest(d_metrics) %>%
          group_by(session_id) %>%
          dplyr::summarize(d_ci_lower = quantile(d, 0.025, na.rm = TRUE),
                           d_ci_upper = quantile(d, 0.975, na.rm = TRUE)) 
      })
    ) %>%
    mutate(d_sig = ifelse((d_ci_lower < 0 & d_ci_upper < 0) | (d_ci_lower > 0 & d_ci_upper > 0), TRUE, FALSE),
           d_ci_width = d_ci_upper - d_ci_lower) %>%
    round_df(3)
  
  data_estimates_with_CIs <- data_D_scores %>%
    full_join(data_D_bootstrapped_CIs, by = "session_id") %>%
    full_join(data_d_scores, by = "session_id") %>%
    full_join(data_d_bootstrapped_CIs, by = "session_id") %>%
    left_join(distinct(select(data_trimmed, session_id, domain), .keep_all = TRUE), by = "session_id") %>%
    select(session_id, domain, 
           D, D_ci_lower, D_ci_upper, D_sig, D_ci_width,
           d, d_ci_lower, d_ci_upper, d_sig, d_ci_width)
  
  # end timer
  end <- Sys.time()
  
  # calculate total time
  end - start

  # save to disk
  write_rds(data_estimates_with_CIs, "models/data_estimates_with_CIs.rds")
  
}

Proportion significant

data_estimates_with_CIs %>%
  arrange(D) %>%
  mutate(ordered_id = row_number()) %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_linerange(aes(x = ordered_id, ymin = D_ci_lower, ymax = D_ci_upper, color = D_sig), 
                 alpha = 0.3) + 
  geom_point(aes(ordered_id, D), size = 0.5) +
  theme_classic() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_color_viridis_d(end = 0.6, direction = -1) +
  xlab("Participants ranked by D score") +
  ylab("IAT D score")

# separated by domains
data_estimates_with_CIs %>%
    group_by(domain) %>%
    arrange(D, .by_group = TRUE) %>%
    mutate(ordered_id = row_number()) %>%
  ungroup() %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_linerange(aes(x = ordered_id, ymin = D_ci_lower, ymax = D_ci_upper, color = D_sig),
                 alpha = 0.3) +
  geom_point(aes(ordered_id, D, color = D_sig), size = 0.5) +
  theme_classic() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_color_viridis_d(end = 0.6, direction = -1) +
  xlab("Participant") +
  ylab("IAT D score") +
  facet_wrap(~domain, ncol = 5)

data_estimates_with_CIs %>%
  summarize(prop_sig = mean(D_sig)) %>% 
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
prop_sig
0.59 
data_estimates_with_CIs %>%
  group_by(domain) %>%
  summarize(prop_sig = mean(D_sig)) %>% 
  arrange(prop_sig) %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
domain prop_sig
David Letterman - Jay Leno 0.32
Lawyers - Politicians 0.37
Artists - Musicians 0.40
Denzel Washington - Tom Cruise 0.41
Gun Control - Gun Rights 0.41
Muslims - Jews 0.41
Meg Ryan - Julia Roberts 0.42
Organized Labor - Management 0.43
Capital Punishment - Imprisonment 0.43
Kobe - Shaq 0.44
Friends - Family 0.44
Mother Teresa - Princess Diana 0.44
Burger King - McDonald’s 0.44
Television - Books 0.44
New York - California 0.46
Receiving - Giving 0.47
Coffee - Tea 0.47
Jazz - Teen Pop 0.48
Microsoft - Apple 0.48
Speed - Accuracy 0.48
Protestants - Catholics 0.49
Private - Public 0.49
Prolife - Prochoice 0.50
Redsox - Yankees 0.50
Wrinkles - Plastic Surgery 0.50
Bill Clinton - Hillary Clinton 0.51
George Bush - John Kerry 0.51
Mountains - Ocean 0.51
State - Church 0.51
Tax Reductions - Social Programs 0.51
Lord of the Rings - Harry Potter 0.52
Pepsi - Coke 0.52
Innocence - Wisdom 0.52
Pants - Skirts 0.52
Team - Individual 0.53
Briefs - Boxers 0.53
Stable - Flexible 0.53
Canadian - American 0.53
Dogs - Cats 0.53
African Americans - European Americans 0.54
Asians - Whites 0.54
Athletic People - Intelligent People 0.54
Helpers - Leaders 0.54
Japan - United States 0.54
Security - Freedom 0.55
Reason - Emotions 0.55
Astrology - Science 0.55
Tall People - Short People 0.55
Protein - Carbohydrates 0.55
Evolution - Creationism 0.56
Corporations - Nonprofits 0.56
Tradition - Progress 0.57
Gay People - Straight People 0.57
Numbers - Letters 0.57
Jocks - Nerds 0.58
West Coast - East Coast 0.58
Single - Married 0.59
Southerners - Northerners 0.59
Hiphop - Classical 0.59
Relaxing - Exercising 0.59
Drinking - Abstaining 0.60
Jews - Christians 0.60
Old People - Young People 0.60
Meat - Vegetables 0.60
50 Cent - Britney Spears 0.61
Rich People - Beautiful People 0.61
Urban - Rural 0.61
Foreign Places - American Places 0.63
Fat People - Thin People 0.63
Traditional Values - Feminism 0.64
Cold - Hot 0.66
Republicans - Democrats 0.67
Atheism - Religion 0.67
Dramas - Comedies 0.68
Realism - Idealism 0.70
Winter - Summer 0.70
Solitude - Companionship 0.71
Technology - Nature 0.72
Career - Family 0.73
Strong - Sensitive 0.73
Effort - Talent 0.76
Conservatives - Liberals 0.77
Determinism - Free will 0.78
National Defense - Education 0.78
Night - Morning 0.79
Past - Future 0.79
Rebellious - Conforming 0.79
Difficult - Simple 0.80
Manufactured - Natural 0.83
Chaos - Order 0.87
Money - Love 0.87
Poor People - Rich People 0.89
Punishment - Forgiveness 0.89
Skeptical - Trusting 0.92
Avoiding - Approaching 0.95 
data_estimates_with_CIs %>%
  group_by(domain) %>%
  summarize(prop_sig = mean(D_sig)) %>% 
  ungroup() %>%
  summarize(min_prop_sig = min(prop_sig),
            max_prop_sig = max(prop_sig),
            mean_prop_sig = mean(prop_sig),
            sd_prop_sig = sd(prop_sig)) %>%
  round_df(2) %>%
  gather() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
key value
min_prop_sig 0.32
max_prop_sig 0.95
mean_prop_sig 0.59
sd_prop_sig 0.14 
trimmed_emprirical_range <- data_estimates_with_CIs %>%
  dplyr::summarize(percentile_025 = quantile(D, 0.025, na.rm = TRUE),
                   percentile_975 = quantile(D, 0.975, na.rm = TRUE)) %>%
  mutate(range_empirical_95_percent = percentile_975 - percentile_025)

41.3% of D scores are not significantly different from zero. So, while they might appear to be relatively large (e.g., D = 0.5), their CI does not exclude zero. Put another way, if we treat the zero point as meaningful, we have insufficient evidence to say whether a given D score represents an IRAP effect in 41.3% of cases in this large sample (11875 participants, 47500 total D scores, 95 domains)).

Half CI width

ggplot(data_estimates_with_CIs, aes(D_ci_width)) +
  geom_density() 

ggplot(data_estimates_with_CIs, aes(D_ci_width, color = domain)) +
  geom_density(adjust = 1.5) +
  theme(legend.position = "none")

library(bayestestR)
point_estimate(data_estimates_with_CIs$D_ci_width/2, centrality = "MAP") %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
MAP
0.38 
data_estimates_with_CIs %>%
  group_by(domain) %>%
  do(point_estimate(.$D_ci_width/2)) %>%
  ungroup() %>%
  round_df(2) %>%
  summarize(min_map_half_ci_width = min(MAP),
            max_map_half_ci_width = max(MAP),
            mean_map_half_ci_width = mean(MAP),
            sd_map_half_ci_width = sd(MAP)) %>%
  round_df(2) %>%
  gather() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
key value
min_map_half_ci_width 0.21
max_map_half_ci_width 0.38
mean_map_half_ci_width 0.36
sd_map_half_ci_width 0.03

Estimation precision

Of course, while treating the zero point as meanining is common, it has been argued to problematic (i.e., by giving rise to false conclusions about the interpretations of the sig, see REF). While the previous analysis represents a useful illustration, it may be more meaningful to consider the precision of estimation of D scores, agnostic to an arbitrary cut-off point.

The width of a D score Confidence intervals was found to be wide: M = 0.66, SD = 0.12, Median = 0.69, MAD = 0.11. Results in the tables above suggest that it doesn’t vary much by trial type or domain. As such, when an individual demonstrates a D score of X, we can more accurately say their D score lies in the range of X ± a median of 0.34.

While the minimum observed D scores was -1.63 and max was 1.75, the outlier scores are clearly visible (see figure XXX). It is therefore likely to be more meaningful to note that 95% of D scores lie within the narrower range of -1.07 to 1.1. It is useful to contextualize the precision of the estimation of a given IAT D score within this total observed ranged of D scores between participants. Specifically, the median CI width noted above represents 31.7% of the (95% trimmed) observed range of D scores. That is to say, the uncertainty around a given D score represents one third of the observed range of D scores: even knowing an individual’s observed D score (e.g., moderately pro-white/anti-black), their ‘true’ D score may lie elsewhere on the range of possible values (e.g., from very pro-white/anti-black to very anti-white/pro-black). An individual IAT sig is therefore quite a poor measure for individual use.

This can also be examined another way by posing the question ‘what proportion of D scores can you tell apart from one another?’ That is, is the probability that a given D score lies outside the CI of all the other D scores’ CIs. For simplicity of implementation, this analysis compares all D scores against all confidence intervals. It is therefore slightly biased by comparing a D score against its own CI as well as all others. However, given the large number of comparisons (i.e., sample size: i.e., 47500 total D scores) this bias is very slight. 95% CIs of this probability value are bootstrapped via case removal and the percentile method using 1000 resamples. The median bootstrapped probability is reported as the estimate for the sake of robustness.

# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_pairwise_comparisons.rds")) {
  
  data_pairwise_comparisons <- read_rds("models/data_pairwise_comparisons.rds")
  
} else {
  
  boots <- data_estimates_with_CIs %>%
    group_by(domain) %>%
    bootstraps(times = 1000)
  
  # helper function to apply workflow to each resample
  helper_function <- function(split) {
    
    estimate <- analysis(split)$D
    ci_lower <- analysis(split)$D_ci_lower
    ci_upper <- analysis(split)$D_ci_upper
    
    n_estimate <- length(estimate)
    n_ci_lower <- length(ci_lower)
    n_ci_upper <- length(ci_upper)
    
    r_estimate <- sum(rank(c(estimate, ci_lower))[1:n_estimate])
    r_ci_upper <- sum(rank(c(ci_upper, estimate))[1:n_ci_upper])
    
    prob_estimate_inferior_to_ci_lower <- 1 - (r_estimate / n_estimate - (n_estimate + 1) / 2) / n_ci_lower
    prob_estimate_superior_to_ci_upper <- 1 - (r_ci_upper / n_ci_upper - (n_ci_upper + 1) / 2) / n_estimate
    
    percent_estimates_inside_cis <- 1 - (prob_estimate_inferior_to_ci_lower + prob_estimate_superior_to_ci_upper)
    
    return(percent_estimates_inside_cis)
    
  }
  
  # apply to each bootstrap
  boot_probabilities <- boots %>% 
    mutate(percent_estimates_inside_cis = furrr::future_map(splits, helper_function)) %>% 
    unnest(percent_estimates_inside_cis) %>%
    select(-splits)
  
  # find CIs using percentile method
  data_pairwise_comparisons <- boot_probabilities %>% 
    summarize(median   = quantile(percent_estimates_inside_cis, 0.500),
              ci_lower = quantile(percent_estimates_inside_cis, 0.025),
              ci_upper = quantile(percent_estimates_inside_cis, 0.975)) %>%
    round_df(3)
  
  
  # save to disk
  write_rds(data_pairwise_comparisons, "models/data_pairwise_comparisons.rds")
  
}

It is not possible to differentiate between two randomly selected D scores in 30.8% (95% CI [30.6, 31]) of cases. This provides additional evidence that IRAP’s individual level precision and therefore clinical utility is low.

Cohen’s d alternative

# plot
data_estimates_with_CIs %>%
  arrange(d) %>%
  mutate(ordered_id = row_number()) %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_linerange(aes(x = ordered_id, ymin = d_ci_lower, ymax = d_ci_upper, color = d_sig), 
                 alpha = 0.3) + 
  geom_point(aes(ordered_id, d), size = 0.5) +
  theme_classic() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_color_viridis_d(end = 0.6, direction = -1) +
  xlab("Participants ranked by Cohen's d score") +
  ylab("IAT Cohen's d score")

data_estimates_with_CIs %>%
  select(session_id, D, d) %>%
  gather(key, value, c(D, d)) %>%
  ggplot(aes(value, color = key)) +
  geom_density()

# separated by domains
data_estimates_with_CIs %>%
    group_by(domain) %>%
    arrange(d, .by_group = TRUE) %>%
    mutate(ordered_id = row_number()) %>%
  ungroup() %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_linerange(aes(x = ordered_id, ymin = d_ci_lower, ymax = d_ci_upper, color = d_sig),
                 alpha = 0.3) +
  geom_point(aes(ordered_id, d, color = d_sig), size = 0.5) +
  theme_classic() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_color_viridis_d(end = 0.6, direction = -1) +
  xlab("Participant") +
  ylab("IAT Cohen's d score") +
  facet_wrap(~domain, ncol = 5)

data_estimates_with_CIs %>%
  summarize(median_half_ci_width = round(median(D_ci_width)/2, 2)) %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
median_half_ci_width
0.34 
ggplot(data_estimates_with_CIs, aes(D_ci_width)) +
  geom_density()

ggplot(data_estimates_with_CIs, aes(d_ci_width)) +
  geom_density()

ggplot(data_estimates_with_CIs, aes(D_ci_width, color = domain)) +
  geom_density() + 
  theme(legend.position = "none")

ggplot(data_estimates_with_CIs, aes(D, D_ci_width)) +
  geom_point(alpha = 0.1) +
  theme_classic()

ggplot(data_estimates_with_CIs, aes(d, d_ci_width)) +
  geom_point(alpha = 0.1) +
  theme_classic()

data_estimates_with_CIs %>%
  group_by(domain) %>%
  summarize(median_half_ci_width = round(median(D_ci_width)/2, 2)) %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
domain median_half_ci_width
50 Cent - Britney Spears 0.34
African Americans - European Americans 0.34
Artists - Musicians 0.36
Asians - Whites 0.35
Astrology - Science 0.35
Atheism - Religion 0.34
Athletic People - Intelligent People 0.36
Avoiding - Approaching 0.26
Bill Clinton - Hillary Clinton 0.35
Briefs - Boxers 0.35
Burger King - McDonald’s 0.36
Canadian - American 0.35
Capital Punishment - Imprisonment 0.35
Career - Family 0.33
Chaos - Order 0.29
Coffee - Tea 0.35
Cold - Hot 0.34
Conservatives - Liberals 0.32
Corporations - Nonprofits 0.35
David Letterman - Jay Leno 0.36
Denzel Washington - Tom Cruise 0.36
Determinism - Free will 0.32
Difficult - Simple 0.30
Dogs - Cats 0.35
Dramas - Comedies 0.33
Drinking - Abstaining 0.35
Effort - Talent 0.33
Evolution - Creationism 0.35
Fat People - Thin People 0.33
Foreign Places - American Places 0.34
Friends - Family 0.36
Gay People - Straight People 0.35
George Bush - John Kerry 0.34
Gun Control - Gun Rights 0.36
Helpers - Leaders 0.35
Hiphop - Classical 0.34
Innocence - Wisdom 0.35
Japan - United States 0.35
Jazz - Teen Pop 0.36
Jews - Christians 0.34
Jocks - Nerds 0.34
Kobe - Shaq 0.35
Lawyers - Politicians 0.37
Lord of the Rings - Harry Potter 0.35
Manufactured - Natural 0.31
Meat - Vegetables 0.33
Meg Ryan - Julia Roberts 0.36
Microsoft - Apple 0.35
Money - Love 0.29
Mother Teresa - Princess Diana 0.36
Mountains - Ocean 0.35
Muslims - Jews 0.36
National Defense - Education 0.31
New York - California 0.36
Night - Morning 0.32
Numbers - Letters 0.35
Old People - Young People 0.34
Organized Labor - Management 0.36
Pants - Skirts 0.35
Past - Future 0.32
Pepsi - Coke 0.34
Poor People - Rich People 0.27
Private - Public 0.36
Prolife - Prochoice 0.36
Protein - Carbohydrates 0.35
Protestants - Catholics 0.36
Punishment - Forgiveness 0.29
Realism - Idealism 0.31
Reason - Emotions 0.36
Rebellious - Conforming 0.32
Receiving - Giving 0.36
Redsox - Yankees 0.35
Relaxing - Exercising 0.34
Republicans - Democrats 0.33
Rich People - Beautiful People 0.35
Security - Freedom 0.35
Single - Married 0.35
Skeptical - Trusting 0.25
Solitude - Companionship 0.33
Southerners - Northerners 0.35
Speed - Accuracy 0.36
Stable - Flexible 0.35
State - Church 0.35
Strong - Sensitive 0.33
Tall People - Short People 0.35
Tax Reductions - Social Programs 0.35
Team - Individual 0.35
Technology - Nature 0.33
Television - Books 0.36
Tradition - Progress 0.35
Traditional Values - Feminism 0.33
Urban - Rural 0.34
West Coast - East Coast 0.35
Winter - Summer 0.34
Wrinkles - Plastic Surgery 0.36 
trimmed_emprirical_range <- data_estimates_with_CIs %>%
  dplyr::summarize(percentile_025 = quantile(D, 0.025, na.rm = TRUE),
                   percentile_975 = quantile(D, 0.975, na.rm = TRUE)) %>%
  mutate(range_empirical_95_percent = percentile_975 - percentile_025)